Stability analysis and numerical evaluations of a COVID-19 model with vaccination

A novel (nonlinear) mathematical model for the transmission of Coronavirus 19 (COVID-19) with eight compartments and considering the impact of vaccination is examined in this manuscript. The qualitative behavior of the system such as the boundedness of solutions, the basic reproduction number, and the stability of the equilibrium points is investigated in detail. Some domestic real data collected from the Kerman University of Medical Science (KUMC) is used to estimate the parameters of the proposed model. We predict the dynamical behavior of the system through numerical simulations based on a combined spectral matrix collocation methodology. In this respect, we first linearize the nonlinear system of equations by the method of quasilinearization (QLM). Hence, the shifted version of Chebyshev polynomials of the second kind (SCPSK) is utilized along with the domain-splitting strategy to acquire the solutions of the system over a long time interval. The uniform convergence and upper bound estimation of the SCPSK bases are proved in a rigorous manner. Moreover, the technique of residual error functions is used to testify the accuracy of the QLM-SCPSK method. The presented numerical results justify the robustness and good accuracy of the QLM-SCPSK technique. The achieved numerical orders of convergence indicate that the QLM-SCSK algorithm has exponential rate of convergence. Using the linearization technique in one hand and the domain-splitting strategy on the other hand, enable us to predict the behaviour of similar disease problems with high accuracy and maximum efficiency on an arbitrary domain of interest.


Introduction
Infectious diseases have always been a limitation of human life.Many epidemics throughout history have killed millions of people.The last of them was the COVID-19 virus, the first case of which was observed in China in December 2019.The disease quickly spread throughout the world within a few months.In the mid of Mach 2020, WHO announced that the infection was an outbreak.A prominent feature of the current pandemic is the high person-to-person transmissibility of the virus, with a basic reproduction number (R 0 ) estimated at 2.2-2.5 in Wuhan, China.This caused an increase in hospitalizations and as a result deaths due to the disease [1,2].Researchers in the field of epidemiology and other branches of biology tried to find treatments or vaccines for this disease.The first vaccines were clinically tested at the beginning of 2021.In the following months, a number of highly effective vaccines entered the market.Mathematical models can describe the process of transmission of infection and its control using vaccination.
Researchers try to understand the dynamics of a disease in the first stage.Hence, they develop control and curing procedures for the diseases.The first modern mathematical models were introduced by Kermak and Mckendrick in 1927 [3].Since 2019, many mathematical models have been presented.Different aspects of the novel COVID-19 have been investigated through mathematical models [4][5][6].Before vaccine production, the effect of quarantine and control factors were investigated [7,8].
Butt et al. developed and analyzed mathematically the extended SEQIHR model.The authors determined the possible control strategies and comprehended the long-term dynamics of disease [9].A SEQIER model introduced to examine the transmission dynamics of COVID-19.The focus is on the effectiveness of hospitalization and quarantine strategies [10,11].A deterministic SEIHR fractional model is developed in [12].The authors performed the mathematical analysis and the design of an optimal control strategy for the proposed Caputo-Fabrizio fractional model.In [13], the authors analyzed a fractional order initial value problem with the Atangana-Baleanu derivative operator to observe transmission dynamics of the infection in a human population.An SVEIR mathematical model was introduced to predict transmission dynamics of COVID-19 infection [14].In [15], the effect of virus rate in the environment in a deterministic model is investigated.Some of the mathematical models of COVID-19 are presented by fractional-order models [16].
Vaccination is very effective to control infectious diseases.Therefore, after producing vaccines for coronavirus, researchers developed dynamical models regarding the effectiveness of vaccines.The authors in [17] developed vaccination strategies for vaccination and investigated the dynamic of their epidemic model.A dynamical model considering the treatment and vaccination saturating function is introduced in [18].The study of the vaccination model under vaccine immunization has been considered in [19].The concerned investigations are devoted to stability theory, numerical simulation, and global-local dynamics [20].One of the items that is investigated in the dynamic analysis of models is bifurcation.Researchers have investigated the existence of bifurcation in different models [21][22][23].
In this paper, a (novel) mathematical model for COVID-19 with an omicron version is presented.In this model, vaccination and its effect on people is considered.In this study, the importance of vaccination to control the epidemic has been investigated.In this model, vaccinated people and those who have not received the vaccine are divided into two groups.Therefore, the rate of infection and hospitalization rate in each group is investigated separately.In this deterministic model, the population is divided into two main parts, people who have been vaccinated and people who have not received any vaccine.It is also assumed that recovered individuals are prepared to reinfection after some time.As follows, some basic dynamic properties of the proposed model are investigated.Using the statistics of the Ministry of Health of Iran in Kerman province, the parameters are estimated over a period of time.
Furthermore, we develop a combined spectral collocation approach based on the shifted version of Chebyshev polynomials of the second kind (SCPSK) to predict the behaviour of COVID-19 disease.To get efficacy, we first transform the nonlinear system of ODEs into a set of linearized system of eight equations to be treated iteratively.Hence, the employed SCPSK basis functions together with the domain decomposition strategy are utilized to find the solutions of the linearized systems.We also prove the convergence of the SCPSK bases and an upper bound estimation for them is performed.It should be noted that spectralbased collocation techniques have been used often to tackle various integral and differential equations due to their simplicity of implementations and high-order accuracy.These spectral methods have been benefited by utilizing numerous basis functions like the Vieta-Lucas [24], Chebyschev [25][26][27], Bessel [28,29], Chelyshkov [30][31][32], and the Schröder [33] polynomials.
Let us illustrate the main achievements of the present work in a concise form as follow: -A novel mathematical model consists of eight equations for the studying of COVID-19 transmission is proposed in which the effect of vaccination is considered.-The parameters of the model are estimated from the real data provided by the KUMC.-The dynamic analysis of the proposed COVID-19 model is performed from theoretical points of stability, boundedness, and the existence of bifurcations.-A combined efficient method called the QLM-SCPSK algorithm based on the technique of quasilinearization along with the Chebyshev spectral collocation approach is designed to solve the given model numerically and confirm the theoretical findings as well.
-The strategy of domain-splitting is further employed to keep the accuracy of the proposed QLM-SCPSK algorithm at a desired level.The exponential convergence of the employed SCPSK bases in the infinity norm is established in a rigourous analysis.
The structure of the current work is provided as follows.A mathematical model considering vaccination is introduced and some basic dynamical properties such as the boundedness of solutions, the derivation of the basic reproduction number, the stability of the equilibrium points, and the existence of bifurcation are investigated in the subsequent "Mathematical model" section.In the next "The Chebyshev polynomials of the second kind: shifted version on [t a , t b ]" section we first give a review of the SCPSK functions.Also, a rigorous mathematical proof is given for the convergence analysis of the SCPSK bases."QLM-SCPSK collocation strategy based on splitting of time interval" section first describes the fundamental ideas of the QLM.Hence, the methodology of the spectral domain decomposition approach is illustrated and finally, the basic steps of the proposed QLM-SCPSK algorithm are given in detail.In "Experimental calculations" section some numerical simulations on a longtime domain are carried out to support the theoretical findings.The last "Conclusions" section includes a summary of the presented research study.

Mathematical model
A dynamical model that is introduced in this paper is based on the Kermak and Mckendrick model [3].The model is introduced for the omicron version of Coronavirus.During the outbreak of the omicron, a large number of people in the community had received the vaccine.Therefore, the population of susceptible, infected, recovered, and isolated are divided into two groups, vaccinated and unvaccinated.Consider, S(t), S v (t) be the populations of unvaccinated and vaccinated susceptible at time t respectively.Denote further by I(t) and I v (t) as the number of unvaccinated and vaccinated infected people.Also, R(t) and R v (t) are unvaccinated recovered and vaccinated recovered persons.Finally, J(t) and J v (t) are the number of unvaccinated and vac- cinated isolated people.
We consider The fraction β of susceptible people who get infected.This fraction of unvaccinated is β ′ .Also, a fraction δ of vaccinated susceptible goes to R v class without infection.The recovery rates are γ 1 and γ 2 for unvaccinated and vaccinated infected people.The unvaccinated and vaccinated infected people are isolated with rates α 1 and α 2 respectively.The recovery rate in J and J v groups are η 1 and η 2 .Two parameters µ 1 and µ 2 are mortality rates because of disease in unvac- cinated and vaccinated isolated populations.In the end, unvaccinated and vaccinated recovered people are reinfected with the rate of θ 1 and θ 2 .In this model, is the birth rate of the population.The dynamical model is given by The flowchart of the model is given in Fig. 1.The above model is subjected to the following initial conditions Let us derive the equilibrium points of the given model Eq.(1).The first equilibrium point is the diseasefree one.It is given by (1) (2) (3) , 0, 0, 0, 0, 0, 0 .The other equilibrium point is obtained as follows where Also, we have defined two parameters x 1 and x 2 as By doing some calculations we find that I * satisfies in the following algebraic equation where the coefficients a i , for i = 0, 1, 2 are given by Also, we have If y 1 > β� , the Eq. 5 has a positive root.Hence, the endemic equilibrium point E * exists.

Proof
Consider the new variable Thus, we have Since, µ < µ i for i = 1, 2 , we arrive at By solving the above linear differential equation in terms of variable X one gets, Now, by tending t → ∞ , we obtain 0 ≤ X ≤ � µ .Therefore, all solutions of the system, i.e., (S, S v , I, I v , R, R v , J , J v ) of the system are confined the region for any ǫ > 0 as t → ∞ .

Basic reproduction number
Let us present the basic reproduction number by R 0 .It is the most important epidemiological parameter.It is defined as the "number of secondary infected individuals caused by a single infected individual in the whole time interval" [34].In this paper, using the method presented in the paper [34], we calculate R 0 value.The compartments that include infected people are I, I v , J , J v .Therefore, we form the following system The above system can be rewritten as where, y = (I, I v , J , J v ) T , and Now, let us obtain the Jacobian matrix related to ϕ and ψ at E 0 as the disease free equilibrium point.We have and As we know that the basic number, R 0 , is the eigen- value of the matrix FV −1 .It is given by some simple calculations as Theorem 2.2 Under the restriction R 0 < 1 , the point E 0 (disease free equilibrium point) is asymptotically stable.

Proof
Denote the Jacobian matrix of COVID-19 system Eq.( 1) by J := J (S, S v , I, I v ) .It has the form where I 1 := I + I v and ζ j := γ j + α j + µ for j = 1, 2 .The former Jacobian matrix at the equilibrium point E 0 can be written as ( 6) .
The characteristic equation is p(x) := det(J (E 0 ) − xI) = 0 .By computing the determinant and after some manipulations we get where we have Owing to the fact that we have R 0 < 1 , one can easily check that all roots are negative.Hence, the equilibrium point E 0 is asymptotically stable.

Bifurcation
In this section, we describe the existence of bifurcation at threshold number R 0 = 1 .We select β to be the bifurcation coefficient that yields bifurcation at R 0 = 1 .

The model takes the form
Here, the transmission rate, β ′ , is less than β .Accord- ing to the available information, consider β ′ = β 30 .Hence, R 0 = 1 gives Next, we construct characteristic polynomial at E 0 and bifurcation parameter β * as follows (8) where q 1 (x) and q 2 (x) are defined in Eq. ( 8).We can eas- ily verify that all roots are negative except one of them is zero.Therefore, there exists bifurcation at R 0 = 1.Now, we compute the right eigenvector w = (w 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 ) T .Corresponding to the zero eigenvalue that satisfies J (E 0 )w = 0.

By solving the equation, we obtain
To compute the left eigenvector Using Theorem 4.1 in [35], we compute the bifurcation coefficients a and b as follow.For the first one, we have where, x 1 = ξ 1 −βS 0 βS 0 v .One can be easily seen that a > 0.
Similarly, for the second parameter b, we have Otherwise, E 0 is unstable.

The Chebyshev polynomials of the second kind: shifted version on [t a , t b ]
In order to devise our multi-domain collocation scheme, we need to introduce the set of (orthogonal) polynomials known as the Chebyshev polynomials of , the second kind (CPSK), see [36].These polynomials on [−1, 1] are defined by the formula where � ∈ [0, π] and j ≥ 0 .The set of CPSK polynomials constitutes an orthogonal system on [−1, 1] with regard to positive function . The shifted version of them are obtained by using p = (2z − 1) on [0, 1].They have the explicit representation as [37] with U 0 (z) = 1 .The second one is obtained as U 1 (z) = 4z − 2 .In addition to Eq. ( 10), we may use the next recursive formula to get the remaining polynomials as Moreover, the orthogonality condition of these shifted polynomials with regard to weight function We finally notice that the roots of U j (z) are located inside the interval (0, 1).Precisely speaking, these zeros are [37] We will utilize the shifted CPSK on an arbitrary interval [t a , t b ] .Toward this end, we take z = (t − t a )/T or t = Tz + t a , where T = t b − t a .So, we have With the preceding transformation, the orthogonality condition Eq. ( 12) is changed accordingly.The orthogonality of the set {U j (t)} ∞ j=0 will be gotten with regard to the weight function w a (t) = T (t − t a ) − (t − t a ) 2 .We have indeed the next formula Finally, based on relations t = Tz + t a and Eq. ( 13) we can locate the zeros of SCPSK U j (t) on (t a , t b ) .The associated zeros are These finite number of points will be utilized as the collocation nodes in the main proposed algorithm, below.

Function approximation: convergence and error analysis of SCPSK
To continue, let us suppose that a given function x(t) is belonged to the weighted L 2 space on [t a , t b ] .It then can be expressed as a summation of SCPSK in the form To obtain the unknown coefficients ν j , j ≥ 0 , the aim is to exploit the orthogonality property Eq.(16).It gives us We first find an upper bound for the coefficients ν j in terms of j.Hence, we show that series form Eq. ( 18) is convergent uniformly on [t a , t b ] .To continue.let assume that M ab := max t∈[t a ,t b ] x ′′ (t) .Moreover, the space L 2 w stands for the space of square-integrable functions on [t a , t b ] with regard to weight function w a (t) defined previ- ously.Thus, we have the next assertion.

Theorem 3.1 For the function
) , which can be written as Eq.(18) we obtain the following estimate

Proof
We begin by the relation Eq. (19) and make the change of variable t − t a = T 2 (1 + cos �) =: r(�) we arrive at (16) We then integrate by parts two times on the last formula.By introducing the auxiliary function we get from Eq. ( 21) the next expression for ν j as We next utilize the fact that x ′′ (r(�)) ≤ M ab and | sin(�)| ≤ 1 .The resulting inequality is A simple calculation can be done to obtain the exact values of the integral on the former line.It follows that By utilizing an odd j > 1 , the resulting integral's value is zero obviously.If we choose a j > 1 to be even, then we have We use the relation j − 1 ≥ j 2 , which holds for all j ≥ 2 .Hence, we reach at the inequality By placing Eq. ( 24) into Eq.( 23) yields the result Eq. (20).
To handle the model problem Eq. ( 1) in practical situations, we require to truncate the infinite series solution Eq. ( 18) into a finite one.If we use only (J + 1) ( J ∈ N ) basis functions, we can approximate x(t) in the way that (21) So, our goal now is to examine the difference between x(t) and x J (t) as an approximation for it.This global error is defined by the formula To get an upper bound for the global error E J (t) , we employ the following relation is valid for the normal CPSK on [−1, 1] given as With the aid of transformation p = 2(t − t a )/T − 1 , we arrive at the same conclusion for U j (t) on [t a , t b ] .Namely, we have The next result establishes that the norm of error �E J � ∞ converges to 0 if we let J goes to infinity.Theorem 3.2 Under assumptions of Theorem 3.1, the error E J (t) converges to 0 as J → ∞ .In fact, we have

Proof
To prove this result, we consider Eq. ( 26) together with Eq. ( 27) to obtain (25) Hence, we apply the upper bound for ν j derived in Eq. (20).The resultant inequality is Utilizing the well-known Integral Test [38] gives us Our proof is established by inserting the foregoing result into Eq.( 28) followed by taking the supremum over all t ∈ [t a , t b ] .By tending J to infinity, we have done the proof.

QLM-SCPSK collocation strategy based on splitting of time interval
Let us emphasize that the spectral matrix collocation approach based on the SCPSK may not yield convergence on a long time interval [t a , t b ] .One remedy is to use a large number of bases on the long domains accordingly to reach the desired level of accuracy.Another approach is to divide the given interval into a sequence of subintervals and employ the proposed collocation scheme on each subinterval consequently.
Towards this end, we split the time interval [t a , t b ] into N ≥ 1 subdomains in the forms Here, we have t 0 := t a and t N := t b .The uniform time step is taken as h = t n+1 − t n = (t b − t a )/N .Note that by selecting N = 1 , we turn back to the traditional spectral collocation method on the whole domain [t a , t b ] .Therefore, on each subinterval K n we take the approximate solution of the model Eq. ( 1) to be in the form Eq. ( 25) as where we utilized the notations as the vector of unknown coefficients and the vector of SCPSK bases respectively.Once we get the all local approximate solutions for n = 0, 1, . . ., N − 1 , the global approximate solution on the given (large) interval [t a , t b ] will be constructed in the form In order to collocate a set of (J + 1) linear equa- tions to be obtained later at some suitable points, we consider the roots of U J +1 (t) on the subinterval K n .By modifying the points given in Eq. ( 17), we take the collocation nodes as At the end, we note that in the proposed splitting approach, the given initial conditions of the underlying model problem are prescribed on the first subinterval K 0 .Once the approximate solution on K 0 = [t 0 , t 1 ] is determined, we utilize it to assign the initial conditions on the next time interval K 1 .To do so, it is sufficient to evaluate the obtained approximation at t 1 .We repeat this idea on the next subintervals in order until we arrive at the last subinterval K N −1 .Below, we illustrate the main steps of our matrix collocation algorithm on an arbitrary subinterval K n for n = 0, 1, . . ., N − 1.

The QLM-SCPSK matrix collocation technique
Our chief aim is to solve the nonlinear COVID-19 system Eq.(1) efficiently by using the spectral method based on SCPSK basis.Towards this end, we first need to get rid of the nonlinearity of the model.This can be done by employing the Bellman's quasilinearization method (QLM) [39].Thus we will get more advantages in terms of running time, especially for large values of J in comparison to the performance of directly applied collocation methods to nonlinear models, see cf. [40][41][42].By combining the idea of QLM and the splitting of the domain we will obtain more gains in terms of accuracy for the approximate solutions of nonlinear model Eq. ( 1).Let us first describe the technique of QLM.For more information, we may refer the readers to the above-mentioned works.
By reformulating the original COVID-19 model Eq. ( 1) in a compact form we get where To begin the QLM process, we assume z 0 (t) is available as an initial rough approximation for the solution z(t) of the COVID-19 system Eq.(31).Through an iterative manner, the QLM procedure reads as follows Here, the notation G z stands for the Jacobian matrix of the COVID-19 system Eq.( 31), which is of size 8 by 8.By performing some calculations we reach the linearized equivalent model form as where M s−1 (t) := J (S s−1 (t), (S v ) s−1 (t), I s−1 (t), (I v ) s−1 (t)) as the Jacobian matrix J previously constructed in Eq. ( 7).

Also we have
Along with the system Eq.( 32) the initial conditions are given due to Eq. ( 2).We now are able to solve the family of linearized initial-value problems Eqs.( 32)-( 33) numerically by our proposed matrix collocation method on an arbitrary (long) domain [t a , t b ] .For this purpose and for clarity of exposition, we restrict our illustrations to a local subinterval K n for n = 0, 1, . . ., N − 1.
In view of Eq. ( 29) by utilizing only ( J + 1 ) SCPSK basis functions, we assume that the eight solutions of system Eq.( 32) can be represented in terms of Eq. ( 29).Thus, we take these solutions at iteration s ≥ 1 as (33) for t ∈ K n .Moreover, by W n,s T we denote the vectors of unknowns for 1 ≤ i ≤ 8 at the iteration s ≥ 1 .Also, the vector of SCPSK basis, i.e., U J (t) is defined in Eq. ( 29).We next provide a decomposition for U J (t) given by Here, the vector Q J (t) including the powers of (t − t n ) introduced by The next object is the matrix F J = (f i,j ) J i,j=0 of size (J + 1) × (J + 1) .The entries of the latter matrix are given in Eq. (15).One can also show that det(F J ) = 0 and it is a triangular matrix.It follows that We then insert the obtained term U J (t) in Eq. ( 35) into Eq.(34).The resulting expansions are (35) We then proceed by nothing that the derivative of the vector Q J (t) can be stated in terms of itself.A vivid cal- culation reveals that From this relation, we are able to derive a matrix forms of the derivatives of the unknown solutions in Eq. (36).
The exact solutions of the linearized system Eq.( 32) can be written in a vectorized form as We next introduce the following block diagonal matrices of dimensions 8(J + 1) × 8(J + 1) as By the aid of the former definitions, the matrix formats of z n J ,s (t) and żn J ,s (t) will rewrite concisely as Here, W n is the successive vector of eight previously defined vector of unknowns ( 37) T .Izadi and Waezizadeh BMC Medical Research Methodology (2024) 24:97 We now can collocate the linearized Eq. ( 32) at the zeros of SCPSK given in Eq. ( 17) on the subdomain K n .We get for s = 1, 2, . . . .Denote the coefficient matrix by M n s−1 and the right-hand-side vector as R n s−1 .These are defined by Let us define further the vectors of unknowns as Consequently, the system of Eq. ( 41) can be stated briefly as and with s = 1, 2, . . . .Before we talk about the fundamen- tal matrix equation, we need to state two vectors Z n s and Żn s in Eq. ( 42) in the matrix representation forms.The proof is easy by just considering the definitions of the involved matrices and vectors in Eq. ( 40).

Lemma 4.1 If two vectors z n
J ,s (t) and żn J ,s (t) in Eq. (40) computed at the collocation points Eq. (30), we arrive at the next matrix forms where the matrix ¯ Q is given by Moreover, two matrices Q, F are defined in Eq. (40).Simi- larly, the vector W n is given in Eq. (40). ( By turning to relation Eq. ( 40) we substitute the derived matrix formats into it.Precisely speaking, after replacing Z n s and Żn s we gain the so-called fundamental matrix equation (FME) of the form where To complete the process of QLM-SCPSK approach, it is necessary to implement the initial conditions in Eq. ( 2) and add them into Eq.(44).So, the next task is to constitute the matrix representation of Eq. ( 2).Let us approach t → 0 in the first relation of Eq. ( 40 This implies that the solution of the model Eq. ( 1) is obtainable on each subdomain K n by iterating n = 0, 1, . . ., N − 1 .On K 0 as the first subdomain, the given initial conditions in Eq. ( 2) will be used to find the corresponding approximations for the system Eq.(1).Hence, this approximate solutions on K 0 evaluated at the starting point of K 1 will be utilized for the initial conditions on K 1 .By repeating this process we acquire all approximations on all K n for 0 ≤ n ≤ N − 1.

REFs and QLM-SCPSK technique
Generally, finding the true solutions of the COVID-19 system Eq.( 1) is not possible practically.In this case, the residual error functions (REFs) help us to measure the quality of approximations obtained by the QLM-SCPSK technique.Once we calculate the eight approximations (44) by the illustrated method, we substitute them into the model system Eq.( 1).In fact, the REFs are defined as the difference between the left-hand side and the right-hand side of the considered equation.On the subdomain K n we set the REFs as for a fixed iteration number s and we have defined L n J ,s := I n J ,s (t) + (I v ) n J ,s (t) for brevity.Analogously, at the fixed iteration s, the numerical order of convergence associated with the obtained REFs can be defined in the infinity norm.These are given by Therefore, the convergence order (Co) for each solution is defined by (46) (47) , ℓ = 1, 2, . . ., 8.

Experimental calculations
We now exploit the proposed QLM-SCPSK collocation technique to solve the COVID-19 system Eq.( 1) numerically.We use Matlab version 2021a installed on a personnel laptop to run our algorithm.Different values of the model parameters will be utilized in the experimental results.Setting s = 5 is sufficient in the iterated QLM to reach the desired accuracy.For each run, to begin computations we use z 0 (t) ≡ 0 as the first rough approximation for the linearized system Eq.(32).
To solve the underlying COVID-19 model, we need to determine the involved parameters in the model.Towards this end, the estimations of these parameters are done in accordance with the statistics given by the Kerman University of Medical Sciences (KUMS) during the period between 22 December 2021 and 5 May 2022.The estimated parameters are listed in Table 1.δ : A fraction of vaccinated susceptible goes to R v class 0.997 The recovery rate of infected people who did not receive the vaccine γ 2 : The recovery rate of infected people who received the vaccine  lines.As one sees the plots obtained by both methods are matched very well.More accuracy is attainable by just increasing the number of bases, namely J, or using a larger number of subintervals N in the computations.Note that for each solution the value of the corresponding population at the endpoint t = 1000 is given for completeness.
As you see, the population of susceptible individuals will decrease in both vaccinated and unvaccinated groups.But with consideration, so a large fraction of the population receives the vaccine, a large number of people directly enter the group of recovered people (R v ) .Hence, the number of infected and isolated individuals is not high compared to the population and a significant decrease is observed in both vaccinated and unvaccinated groups.According to the prediction of these graphs, with Fig. 2 The behavior of susceptible population: Unvaccinated S(t) (left) and vaccinated S v (t) (right) obtained via QLM-SCPSK methodology using J = 10 and N = 100 on [0, 1000] Fig. 3 The behavior of infected population: Unvaccinated I(t) (left) and vaccinated I v (t) (right) obtained via QLM-SCPSK methodology using J = 10 and N = 100 on [0, 1000] the passage of time, the population goes through an epidemic, which is consistent with the reality that happened in the province of Kerman.
We next show the effect of utilizing a diverse number of bases, J = 10, 20 , on the computed solutions.In this respect, we compute the REFs formulae in Eq. ( 46) related to the approximate solutions of the COVID-19 model Eq.(1).To save space, we only visualize the REFs corresponding to the susceptible and recovered populations.These results for both vaccinated and unvaccinated counterparts are displayed in Figs. 6 and 7.Not that we used the fixed number of subintervals as N = 100 .It can clearly seen that the desirable form of accuracy is attainable if one increases the number of bases accordingly.
As we mentioned, we can also get the smaller magnitude of REFs by increasing the number of subintervals N in the proposed domain decomposition spectral QLM-SCPSK approach.Towards this end, let use take N = 50, 100 and N = 200 .For a fixed J = 15 , we show the results of REFs associated with the recovered populations (vaccinated and unvaccinated) in Fig. 8. From the plots shown in Fig. 8, it can be obviously concluded that high-order accuracy is attained for smaller lengths of intervals.
Finally, we compute L ∞ ℓ error norms and the related numerical order of accuracy as shown by Co ℓ J for ℓ = 5, 6, 7, 8 , namely for recovered and isolated popula- tions.The results by using diverse values of J = 2 k , for Fig. 4 The behavior of recovered population: Unvaccinated R(t) (left) and vaccinated R v (t) (right) obtained via QLM-SCPSK methodology using J = 10 and N = 100 on [0, 1000] Fig. 5 The behavior of isolated population: Unvaccinated J(t) (left) and vaccinated J v (t) (right) obtained via QLM-SCPSK methodology using J = 10 and N = 100 on [0, 1000] Fig. 6 The REFs related to susceptible population: Unvaccinated R n 1,J (t) (left) and vaccinated R n 2,J (t) (right) obtained via QLM-SCPSK methodology using J = 10, 20 and fixed N = 100 on [0, 1000] Fig. 7 The REFs related to recovered population: Unvaccinated R n 5,J (t) (left) and vaccinated R n 6,J (t) (right) obtained via QLM-SCPSK methodology using J = 10, 20 and fixed N = 100 on [0, 1000] Fig. 8 The REFs related to recovered population: Unvaccinated R n 5,J (t) (left) and vaccinated R n 6,J (t) (right) obtained via QLM-SCPSK methodology using a fixed J = 15 and various N = 50, 100, 200 on [0, 1000] Izadi and Waezizadeh BMC Medical Research Methodology (2024) 24:97 k = 1, 2, . . ., 5 are listed in Table 2.For this purpose, we use N = 10 and the computational interval is set as [0, 100].By looking at the results tabulated in this table one can infer that the proposed spectral QLM-SCPSK approach has an exponential order of accuracy.

Parameter study
In this part, the relationship between β, γ and R 0 is exam- ined.As you can see in the next Fig. 9 (left picture), as the γ value increases, R 0 decreases.The lowest value of R 0 is at γ = 0.3 .Therefore, if the recovery rate reaches 0.3, the epi- demic will disappear.The relationship between β and R 0 is shown in the same figure but on the right panel.As you see, there is a direct relationship between them.

Conclusions
A novel mathematical model for studying the COVID-19 pandemic disease has been suggested by dividing the population into vaccine and nonvaccine groups.From a dynamical point of view, the boundedness of the system was proved and the basic reproduction number was obtained to prove the stability of the equilibrium points.Moreover, the existence of bifurcation for the COVID-19 system was discussed.To get the solutions of this system, an efficient spectral method based on the shifted Chebyshev polynomials of the second kind (SCPSK) combined with the quasilinearization methodology (QLM) was used.The methodology of domainsplitting further was implemented to keep the accuracy of the proposed QLM-SCPSK approach at a desired level on a long time interval.The convergence analysis of the SCPSK basis functions in the L ∞ norm was per- formed and an upper bound estimation for the error was done.The involved parameters were estimated by the real data provided by the Kerman University of Medical Sciences during the period between 22 December 2021 and 5 May 2022.Numerous computational experiments based on the proposed QLM-SCPSK technique were conducted to predict the behavior of disease over a long time interval.The presented QLM-SCPSK technique with high accuracy and efficacy can be generalized to solve similar epidemiological models with integer-order and fractional-order derivatives and even with more than eight equations.

Fig. 1
Fig. 1 Schematic view of the COVID model

Definition 3 . 1
The shifted CPSK (SCPSK) on [t a , t b ] is defined via relation From the above change of variable, we are able to represent the summation Eq. (10) as where(9) U j (p) := sin [(j + 1)�]/ sin �, p = cos �,(10) ).It gives us We then replace eight rows of the augmented matrix [B n ; R n s−1 ] by the already obtained row matrix [B 0,n ; R n s−1,0 ] .Denote the modified FME by The initial conditions are taken asTo start computations, we first set J = 10 and the com- putational domain is taken as [t a , t b ] = [0, 1000] in terms of days spent during the pandemic.The number of subintervals used is N = 100 .Through Figs. 2, 3, 4 and 5 we show the behaviors of the approximate solutions obtained via QLM-SCPSK procedure for susceptible, infected, recovered, and isolated populations in both unvaccinated and vaccinated parts.These curves are presented by dashed red lines.Besides, the outcomes of the well-known Matlab function ode45 are also depicted to validate our approximations.These plots are shown by solid black (48) S 0 = 1184950, S v0 = 1815050, I 0 = 760, I v0 = 25, R 0 = 2720, R v0 = 25, J 0 = 152, J v0 = 5.

Table 1
Parameters: Descriptions and their values

Table 2
Maximum value of REFs and the associated Co ℓ J for ℓ = 5, 6, 7, 8 obtained by using the QLM-SCPSK approach with N = 10 and various J on [0, 100]